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Continuous-time determinantal algorithm is proposed for the quantum Monte Carlo simulation 
of the interacting fermions. The scheme does not invoke Hubbard-Stratonovich transformation. 
The fermionic action is divided into two parts. One of them contains the interaction and certain 
additional terms; another one is purely Gaussian. The first part is considered as a perturbation. 
Terms of the series expansion for the partition function are generated in a random walk process. 
The sign problem and the complexity of the algorithm are analyzed. We argue that the scheme 
should be useful particularly for the systems with non-local interaction. 

The family of quantum Monte Carlo (QMC) methods is the most universal tool for the numerical study of quantum 
many-body systems with strong correlations. So-called determinental methods 0, Q are commonly used For the 
interacting fermions, because other known QMC schemes (like Stochastic Series Expansion or worm algorithms 
0) suffer an unacceptably bad sign problem for this case. The partition function of the system is presented as a 
sum over an ensemble of the systems of uncoupled fermions. This idea was proposed by Scalapino and Sugar pj. 
Their algorithm was later significantly improved by Hirsh Q • The scheme by Hirsh is now almost standard for QMC 
simulations. Two points are essential for the method: first, the imaginary time is artificially discretized, and then 
a discrete Hubbard-Stratonovich transformation is performed at each time slice to decouple the fermionic degrees of 
freedom. After this decoupling, the fermionic degrees of freedom can be integrated out, and Monte Carlo sampling 
should be performed in the space of auxiliary Hubbard-Stratonovich spins. For a system of N atoms the number of 
spins scales as (3N/tq for the local (short-range) interaction and as (3N 2 /r for the long-range one (f3 is the inverse 
temperature, and r is a time slice). 

One can point out the following drawbacks of the Hirsh algorithm. First, for large systems and small temperature, 
the sign problem results in the exponential growth of computational time. This is a principal problem of QMC, and 
there is no solution on the horizon now. However, for relatively small clusters, and in particular for dynamical mean- 
field calculations Q, the sign problem is not of crucial importance. Further, the non-locality of interaction hampers 
the calculation, because it is hard to simulate systems with a large number of auxiliary spins. It is nearly impossible to 
simulate the system with interactions that are non-local in time, when the number of spins is proportional to ((3/tq) 2 . 
And finally, the time discretization leads in a systematic error of the result. 

In 1999th a continuous-time modification of the algorithm was proposed 0. It was based on a series expansion for 
the partition function in powers of interaction. The scheme is free of systematic errors. The Hubbard-Stratonovich 
transformation is still invoked. The number of auxiliary spins scales similarly to the discrete scheme. 

In this communication we present a continuous-time algorithm which is based on the series expansion as well but 
does not deal with any auxiliary-field variables. For the case of Hubbard model the sign problem and the computational 
time are found to be similar to what occurs for the Hirsh scheme. The principal advantages of the present algorithm 
are related to the different scaling of the computational time for non-local interactions. The paper is aimed at a 
general description of the algorithm and the estimation of the computation time. Our pilot numerical results are not 
presented here because the systematic study is not yet done. 

Consider a system of fermions with pair interaction. For the most general situation, the partition function can be 
presented as follows: 

Z = SpTe- s (1) 
S = JJ t r r 'cl,c r drdr' +JJJJ w r ^\c\,c^c\,c^dr 1 dr l 1 dr 2 dr' 2 . 

Here T is a time-ordering operator, r — {r, s,i} is a combination of the continuous imaginary-time variable r, spin 
orientation s and the discrete index i numbering the single-particle states in a lattice. For the single-band models i is 
an atom number. For electrons s =| or j. Integration over dr implies the integral over r, and the sum over all lattice 
states and spin projections: J dr = dr. 

We borrow the linear-algebra style for sub- and superscripts to make the notation clearer. The creation and 
annihilation operators are labelled as covariant and contravariant vectors, respectively. The labelling for coefficients 
t, w is chosen to present all integrands like scalar products of tensors. In principle, we could declare the summation 
over the repeating indices and would not explicitly write integrals. However, we keep them for convenience. 
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Let us introduce an additional quantity a r r ,, and split S into two parts. Up to an additive constant 

S = So + W, (2) 
So =JJ (V +// a r r l(w r rX 2 +«4r Jdradr'a) 4,c r drdr', 

W = JJJJ Wr\^{c\,c^ -afyicl^ - a^)dridr , 1 dr 2 dr , 3 . 

The freedom to choice a r r , is to be used later to optimize the algorithm, particularly to minimize the sign problem. 

We consider So as an unperturbed action and switch to the interaction representation. The perturbation-series 
expansion for Z is as follows: 

z = E feLo / dr i / dr i ■ ■ ■ I dr 2k J dr' 2k fl k (n , r[ , . . . , r 2k , r' 2k ) (3) 

O, - 7, (- 1 )* , „. r i r 2 . „,, r 2k-l r 2k n r ir2 ...r 2k 

' 'l'2-'-'2fc 

where Zo is a partition function for the unperturbed system and 

D r r r~ r r ?" =< T(c\,c^ - a r r j) ■ ... • (4, c r » - a r r l») > . (4) 

2k ' 1 ' 1 ' 2fc ' 2fc 

Hereafter the triangle brackets denote the average over the unperturbed system, < A >= Z _1 Sp (TAe~ s °). 

To obtain an expression for D, one can start from the case of a = 0. Since So ls Gaussian, Wick theorem applies. 
Therefore D is a determinant 2k x 2k matrix. Two-point Green functions 

g r r , =< Tc\,c r > (5) 

form this matrix for a = 0. Now, the non-zero a can be taken into account. It is easy to prove that 

DXr-- r r ? k =det\\g r r i-S(i-j)a r ;,\\ (6) 

1 2" - 2k j j 

where 8(i — j) is a delta-symbol. 

Like any other other QMC scheme, the proposed one is based on a Markov process. Points in the configuration 
space are determined by the number k and the set ri,r[, r 2/ t, r' 2k . Suppose for a moment that Q is always positive, 
and consider a random walk with a probability of Z~ 1 £lk{r\, r' 1; r 2 k, f^k) t° visit each point. Then for example the 
Green function can be calculated from 

Z- 1 < Tcl,c r e- W >= ff ;,(r 1 ,r' 1 ,...,r^) = (7) 
= Efc / dr i I dr'i-J dr 2k gr<(ri,r[, ...,r^.)fi fe (ri,ri, -,r' 2k ), 
where <?£'( r i> r i> ■■■' r 2fc) determines the Green function for a current realization 

ffPCri.ri, ...y 2k ) = (D^ll)- 1 < Tcl,c r (cl,c^ - aQ) ■ ... • {c\,c^ -<?)>. (8) 

The overline stands for the averaging over the above-mentioned random walk. The important notice is that the series 
expansion for an exponent always converges, therefore the discussed averaging is always well-defined in a mathematical 
sense. 

The sign problem should be discussed first. If f2 is not always positive, one should work with a probability |f2| 
and calculate <?sign(f2)/sign(f2) instead of g. For a practical calculation it is desirable to maximize the average sign, 
because otherwise a computational error bar is unacceptable. 

A proper choice of a can suppress the sign problem in certain cases. To be concrete, let us consider a Hubbard 
model. In this model the interaction is local in time and space, and only electrons with opposite spins interact. 
Therefore it is reasonable to take ct^],^ — S(t — r')5(i — i') a T> the similar for a^, and aj = aj = 0. The perturbation 
W becomes 

IV r Hubbard = U /(n-f(r) - af)(nj.(r) - ajdi (9) 

Here Hubbard U and the occupation number operator n — c^c are introduced. The Gaussian part of Hubbard action 
is spin-independent and does not rotate spins. This means that only g\,g\ do not vanish, and the determinant in © 
is factorized 



D r rr~ r ? k = DTf-l k - 1 Dl r t-l k = £> T £>i (10) 

r, rU...rL, r i r 3~- r 2k — 1 r ^ r -> ■■■ r ^^ 1 
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For the case of attraction U < one should choose 

a T = "1 = a > ( n ) 

where a is a real number. For this choice gj = gj, and therefore £)j = D^. f2 is always positive in this case, as it 
follows from formula 

This choice of a is useless for a system with repulsion, however. Compared to the case of attraction, another sign 
of w at a-f = aj, results in the alternating signs of flk with odd and even k (£(• Another condition for a is required. 
The particle-hole symmetry can be exploited for the Hubbard model at half-filling. In this case, the transformation 
ct — ► cj, converses the Hamiltonian with repulsion to the same but with attraction. Therefore the series © in powers 
of W = U J(ti^(t) — a)(hi(r) — a)dt does not contain negative numbers, in accordance to the previous paragraph. 
The value of trace in Jjyl is independent of a particular representation. In the original (untransformed) basis the above 
written W reads as U J (n^(r) — a)(n|(r) — 1 + a)dt. We conclude that 

a-f = 1 — cx-i = a 

eliminates the sign problem for repulsive systems with a particle-hole symmetry. 

Of course, the average sign for a system with repulsion is not equal to unity in a general case. As usual the 
exponential fall-off occurs for the large systems or small temperature. 

Now we discuss how to organize a random walk in practice. We need to perform a random walk in the space of 
k; ri, r^, r2k,r' 2k . Two kinds of trial steps are necessary: one should try either to increase or to decrease k by 1, and, 
respectively, to add or to remove the four corresponding "coordinates". A proposition for r^k+ii r 2fc+i' r 2k+2, r 2k+2 
should be generated for the "incremental" step. The normalized modulus 

iMr>&£&: a a i (13) 

\\w\\ = ffff\w r r 'x'\drdRdr'dR' 

can be used as a probability density for this proposition. Then the standard Metropolis acceptance criterion can be 
constructed using the ratio 



HI 
fc + i 



n— r 2 k+2 

(14) 



D') 

r l"- r 2fc + 2 



For the " decremental" step one should take a random integer j between and k — 1 and consider a probability of 
the removing r2j+i , r 2j+i > r 2j+2 , r 2j+2 ■ The detailed balance condition should be satisfied, therefore an inverse of (|14|) 
should be used in the acceptance criterion. 

It is possible to introduce more complicated trial steps in a similar way. It can be an increment or decrement of k 
by 2, shifts of r etc. Normally, the basic steps discussed above should be sufficient. We shall demonstrate however 
that for certain cases the increments (decrements) by 2 are required. 

The most time consuming operation of the algorithm is a calculation of the ratio of determinants. This task is 
relatively easy if Gaussian part of action is described by a Hamiltonian So — J drHo (r) and the interaction is also a 
combination of the local in time parts: 

Wrlrl OC <5(Ti - t[)S(t 2 - Tq). (15) 

One can call this system " almost Hamiltonian" . Purely Hamiltonian systems and particularly Hubbard model also 
belong to this class, of course. For this case the expression Q for D can be written explicitly in the Schrodinger 
representation: 



DJ 



ri..-r 2 k 
■ ' ...r' . 



Z^Sp (e- TlHo (c{ c 11 - a i , 1 )e (Tl ~ T2)ffo • ... • (cj ,c i2fc - aj k )e T2hH °) , (16) 

\ 11 2fc 2k / 



where T\...T2k are pre-ordered in time. Further, an identity c\,c l — a = — aexp(^c|,c i ) holds with £ = — a -1 for i ^ i' 
and ^ = ln(l — a -1 ) for i = i'. Trace here can be calculated in a Scalapino-Sugar manner 0, that requires the 
calculation of oc N x N matrix determinant for a system of N atoms. Although the calculation of a determinant from 
scratch takes oc N 3 operations, the fast-update technique can be used here, resulting in a oc N 2 operation-count. 
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For the general situation, formula © is to be used to estimate D. Here we present an estimation for the typical 
value of fc. An average value of l|14(l determines an acceptance rate for QMC sampling. It is reasonable to expect that 
by the order of magnitude this rate is not much less then unity, and consequently 



IM 



(17) 



D') 

T l---' 2k + 2 



D r J 



If the sign problem is absent and the sign of w r r \rl is the same for all arguments, the right-hand side can be interpreted 
as an expectation value of \W\ (compare with formula 0), i.e. 

k^]W\. (18) 

For the Hubbard lattice of N atoms, for instance, \W\ oc j3\U\N. Using fast updates one can make a step in oc fc 2 
operations. 

To obtain the fastest procedure one needs to minimize fc with respect to a's. It should decrease not only the 
computational time for a single random-move, but the autocorrelation length of the random walk as well. Therefore 
the minimization of k is desirable for the almost-Hamiltonian systems also. 

For example, once the the condition (|12|l for the Hubbard model with attraction is fulfilled, the value of a still can 
be adjusted. We obtain 



\W\ oc n T (r)n x (r) - 2na + a, (19) 

where n is the average filling factor. Minimization gives a physically reasonable suggestion a = H. 

It is useful to analyze a toy single-atom Hubbard model to get a filling of the behavior of series © . The two parts 
of action are 

So = J(- f i + C/a x )n T (r) + (-// + C/a T )n x (r))dr; (20) 
w = U J(ni(r) - ai)(n T (T) - a T )dr. 

Here [i is a chemical potential. For a half-filled system \i = U/2. For this model, it is easy to calculate traces in JTBJ, 
and 10 becomes 



{-Uo^aCf ( 1 + e fl(M-t/ai) ( i_ a -i ) fc) (l + e/^-t^i-a- 1 )*) (21) 



Consider a case of repulsion (U > 0). Let us use a condition (|12|l for an arbitrary filling factor. The later expression 
can be presented in the form 



nk = e 0(n-u a )(U?_)_ (i + e /J(M-tf+t>a)(i_ a -i)*) (l + e^+^a-a- 1 )*) (22) 

For fj, = U/2 the value of fi/t is positive for any a. For a general filling factor, the situation depends on the value of 
a. For < a < 1 negative numbers can occur at certain fc. Outside this interval all terms are positive, and there is 
no sign problem for a single-atom system under consideration. 

Since the sign problem persists already for the impurity problem for < a < 1, such a choice is also not suitable 
for the TV-atom repulsive Hubbard system. On the other hand, minimization of W requires a to be as close to this 
interval as possible. Therefore it is reasonable to take a = 1 or slightly above. This is the same as zero or a small 
negative value, since a-f = 1 — a^. In this case our pilot calculations and the discrete-time scheme demonstrate 
approximately the same sign problem. 

There is also an important note about the system at half-filling. For this case the optimal choice is a-\ = ay = 1/2. 
In this case one can observe from (1221 that Slfc = for any odd fc. This is a consequence of a particle-hall symmetry. 
This is an example of the situation when the trial increase/decrese of fc by 1 does not help, and it is necessary to 
consider changes of fc by 2. 

In a conclusive part, we would like to discuss possible benefits of the proposed algorithm over another determinantal 
schemes. We demonstrated that for a Hubbard-type models the computational time for a single trial step scales as 
(PUN) 2 for the general case and as iV 2 for a Hamiltonian system. This is the same scaling as for the schemes based 
on a Stratonovich transformation. An important difference occurs however for the non-local interactions. Consider, 
for example, a system with a large Hubbard U and much smaller but still important Coulomb interatomic interaction. 
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One needs to introduce iV 2 auxiliary fields per time slice instead of N to take the long-range forces into account. This 
slows the calculation dramatically. On the other hand, the complexity of the present algorithm should remain almost 
the same as for the local interactions, because \W\ does not change much. It should be also possible to study the 
interactions retarded in time, particularly the effects related to dissipation. 

The proposed algorithm operates in a continuous time and does not invoke Statonovich transformation. A method- 
ologist can ask if there is a similar scheme for a discrete time. The answer is yes, one should only replace integrals 
with sums in all formulae. The time-scaling should remain the same as for the continuous time. The discrete-time 
formulation might be handier for the practical realization. 

I am grateful to A.I. Lichtenstein and F. Assaad for their very valuable comments. This research was supported in 
part by the National Science Foundation under Grant No. PHY99-07949. 
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